DEMO: H wavefunctions#
Requirements | Importing libraries#
import matplotlib, matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
from scipy.constants import physical_constants
from matplotlib import cm, colors
import plotly.graph_objects as go
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
try:
from google.colab import output
output.enable_custom_widget_manager()
print('All good to go')
except:
print('Okay we are not in Colab just proceed as if nothing happened')
Okay we are not in Colab just proceed as if nothing happened
1. Describing a radial function \(R_{nl}(r)\)#
\[R_{nl} = \sqrt{ \Big( \frac{2}{n a_0} \Big)^3 \frac{(n-l-1)!}{2n(n+l)!}} \cdot e^{-\frac{r}{n a_0}} \cdot \Big( \frac{2 r}{n a_0} \Big)^l \cdot L^{2l+1}_{n-l-1}\Big( \frac{2 r}{n a_0} \Big)\]
def radial_function(r,
n=1,
l=0):
'''Rnl(r) normalized radial function
r: radial distance Float (0, inf)
n: principal quantum number Int (1,2,3... inf)
l: angular quantum number Int (0,1,2,... n-1)
'''
from scipy.special import genlaguerre
a0 = 1 # Bohr radius equal to 5.29e-11 m
prefactor = np.sqrt( ((2 / n * a0) ** 3 * (np.math.factorial(n - l - 1))) / (2 * n * (np.math.factorial(n + l))) )
laguerre = genlaguerre(n - l - 1, 2 * l + 1)
p = 2 * r / (n * a0)
return prefactor * np.exp(-p / 2) * (p ** l) * laguerre(p)
r = np.linspace(0, 20, 1000)
Rnl = radial_function(r, n=2, l=0)
#plt.plot(r, Rnl)
#plt.plot(r, Rnl**2)
plt.plot(r, r**2 * Rnl**2)
plt.xlabel(r'$r [a_0]$')
plt.ylabel(r'$R_nl(r) r^2$')
/tmp/ipykernel_15473/1039204017.py:14: DeprecationWarning: `np.math` is a deprecated alias for the standard library `math` module (Deprecated Numpy 1.25). Replace usages of `np.math` with `math`
prefactor = np.sqrt( ((2 / n * a0) ** 3 * (np.math.factorial(n - l - 1))) / (2 * n * (np.math.factorial(n + l))) )
Text(0, 0.5, '$R_nl(r) r^2$')
2. Describing an angular function | Spherical harmonic \(Y_{l}^{m}(\theta, \varphi)\)#
We will make use of sph_harm function from Scipy
We can also build up spherical harmonics using Associated Legendre function using lpvm
\[
Y_{lm}(\phi,\theta) = \sqrt{\frac{2l+1}{4\pi} \frac{(l-m)!}{(l+m)!} } P_{lm}(cos \phi) \cdot e^{im\theta}
\]
from scipy.special import sph_harm
sph_harm(0, 0, np.pi, np.pi) # test a few spherical harmonics
(0.28209479177387814+0j)
def plot_Yml(l, m):
'''Visualizing spherical harmonics using sph_harm funcion from Scipy special function library
Note that the name of angles is different from the notation adopted in QM textbooks!
l: angular quantum number (0,1,2,...)
m: magnetic quantum number (-l,...l)
'''
#Creade grid of phi and theta angles for ploting surface mesh
phi, theta = np.linspace(0, np.pi, 100), np.linspace(0, 2*np.pi, 100)
phi, theta = np.meshgrid(phi, theta)
#Calcualte spherical harmonic with given m and l
Ylm = sph_harm(m, l, theta, phi)
R=abs(Ylm)
# Let's normalize color scale
fcolors = Ylm.real
fmax, fmin = fcolors.max(), fcolors.min()
fcolors = (fcolors - fmin)/(fmax - fmin)
# Since we want to plot on cartesian reference frame we will use cartesian coordiniates x, y, z using R as the absolute value of Yml
# Try a plot without R part.
fig = go.Figure(data=[go.Surface(x=R*np.sin(phi) * np.cos(theta),
y=R*np.sin(phi) * np.sin(theta),
z=R*np.cos(phi),
surfacecolor=fcolors,
colorscale='balance')])
# Show the plot
fig.update_layout(title=fr'$Y_{l, m}$', autosize=False,
width=700, height=700,
margin=dict(l=65, r=50, b=65, t=90)
)
fig.show()
plot_Yml(l=1, m=0)
3. Describing the normalized probability as wavefunction squared \(|\psi _{nlm}(r,\theta ,\varphi)|^2\)#
\[\psi_{nlm} = R_{nl}(r) \cdot Y_{l}^{m}(\theta, \varphi)\]
def normalized_wavefunction(n,
l,
m,
max_r = 10):
'''Ψnlm(r,θ,φ) normalized wavefunction
by definition of quantum numbers n, l, m and a bohr radius augmentation coefficient'''
# Set coordinates grid to assign a certain probability to each point (x, y) in the plane
x = y = np.linspace(-max_r, max_r, 501)
x, y = np.meshgrid(x, y)
r = np.sqrt((x ** 2 + y ** 2))
# Ψnlm(r,θ,φ) = Rnl(r).Ylm(θ,φ)
psi = radial_function(r, n, l) * sph_harm(m, l, 0, np.arctan(x / (y + 1e-10)))
return np.abs(psi) ** 2
4. Plotting wavefunction electron probability density plots#
Light shaded areas in the orbital cross-sections represent a high probability of a particle being present in that region.
def plot_orbital_2d(n,l,m, max_r = 10):
fig, ax = plt.subplots()
im = ax.contour(normalized_wavefunction(n, l, m, max_r=max_r), 16, cmap='rocket', extent=[-max_r, max_r,-max_r, max_r])
ax.set_title(r'$|\psi_{{({0}, {1}, {2})}}|^2$'.format(n, l, m), fontsize=15)
ax.set_aspect('equal')
plot_orbital_2d(5,3,0, max_r=20)
/tmp/ipykernel_15473/1039204017.py:14: DeprecationWarning:
`np.math` is a deprecated alias for the standard library `math` module (Deprecated Numpy 1.25). Replace usages of `np.math` with `math`